home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / testmatrix / cauchy.r < prev    next >
Text File  |  1994-12-20  |  2KB  |  63 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Cauchy matrix.
  4.  
  5. // Syntax:    C = cauchy ( X , Y )
  6.  
  7. // Description:
  8.  
  9. //    Compute the Cauchy matrix, where X, Y are N-vectors, is the
  10. //    N-by-N matrix with C[i;j] = 1/(X[i]+Y[j]).  By default, Y =
  11. //    X. Special case: if X is a scalar CAUCHY(X) is the same as
  12. //    CAUCHY(1:X). 
  13.  
  14. //    Explicit formulas are known for DET(C) (which is nonzero if X
  15. //    and Y both have distinct elements) and the elements of
  16. //    INV(C). C is totally positive if 
  17. //        0 < X(1) < ... < X(N) and 0 < Y(1) < ... < Y(N).
  18.  
  19. //      References:
  20. //      D.E. Knuth, The Art of Computer Programming, Volume 1,
  21. //           Fundamental Algorithms, second edition, Addison-Wesley, Reading,
  22. //           Massachusetts, 1973, p. 36.
  23. //      E.E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications,
  24. //           Linear Algebra and Appl., 149 (1991), pp. 1-18.
  25. //      O. Taussky and M. Marcus, Eigenvalues of finite matrices, in
  26. //           Survey of Numerical Analysis, J. Todd, ed., McGraw-Hill, New York,
  27. //           pp. 279-313, 1962. (States the totally positive property on p. 295.)
  28.  
  29. //    This file is a translation of cauchy.m from version 2.0 of
  30. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  31. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  32.  
  33. //-------------------------------------------------------------------//
  34.  
  35. cauchy = function ( x , y )
  36. {
  37.   local (x, y)
  38.  
  39.   n = max(size(x));
  40.  
  41.   //  Handle scalar x.
  42.  
  43.   if (n == 1)
  44.   {
  45.     n = x;
  46.     x = 1:n;
  47.   }
  48.  
  49.   if (!exist (y)) { y = x; }
  50.  
  51.   x = x[:]; 
  52.   y = y[:];   // Ensure x and y are column vectors.
  53.  
  54.   if (any(size(x) != size(y))) {
  55.     error ("Parameter vectors must be of same dimension.");
  56.   }
  57.  
  58.   C = x*ones(1,n) + ones(n,1)*y.';
  59.   C = ones(n,n) ./ C;
  60.  
  61.   return C;
  62. };
  63.